library(tidyverse)


Municipality_data <- read.csv("Data/Municipality_data_2021_01_18.csv") %>%
  as_tibble() %>%
  filter(AnnéeDécès != "") %>%
  # Filtering what is actually population data
  filter(AnnéeDécès > 1900) # Filtering out 3954907 MILAN RENEE/     2           1923             6             7 99                              127 SAINT-EUF┬ÉENIA      …        201
Population_data_extrapolated <- readRDS("Results/Population_data_extrapolated.rds")


selected_years <- c(2015:2017, 2019, 2020)

mortality_period_municipality_level <- function(start_date, end_date) {
  Municipality_data %>%
    filter(as.numeric(JourDécès) != 0) %>%
    # Excluding 15238 death records
    # Filtering what is actually population data
    mutate(Death_date = lubridate::ymd(paste(AnnéeDécès, MoisDécès, JourDécès, sep = "/"))) %>%
    filter(Death_date >= start_date) %>%
    filter(Death_date <= end_date) %>%
    mutate(Municipality_code = Code) %>%
    group_by(Municipality_code) %>%
    summarise(NbDeath = n())
}

Mortality_data <- selected_years %>%
  lapply(function(y) {
    mortality_period_municipality_level(start_date = ymd(paste(y, "/03/15", sep = "")), end_date = ymd(paste(y, "/04/30", sep = ""))) %>%
      mutate(year = y)
  }) %>%
  bind_rows()

Mortality_data_joined <- Population_data_extrapolated %>%
  select(Code, Population, year) %>%
  filter(year %in% selected_years) %>%
  rename(Municipality_code = Code) %>%
  left_join(Mortality_data) %>%
  mutate(NbDeath = ifelse(test = is.na(NbDeath), yes = 0, no = NbDeath))

selected_cities <- Mortality_data_joined %>%
  filter(year == 2020) %>%
  filter(Population >= 710) %>%
  # Filtering which should match the electoral data filtering
  .$Municipality_code

full_data <- Mortality_data_joined %>%
  dplyr::filter(Municipality_code %in% selected_cities) %>%
  mutate(municipality_idx = Municipality_code %>% as.factor() %>% as.numeric())

res <- readRDS("Results/vb_estimated_shrunk_hazard_rate_2021_01_18_rep2.rds")
vb_fit <- res$vb_fit
fit_data <- res$data

MAP_fit <- readRDS("Results/MAP_estimated_shrunk_hazard_rate_2021_01_18_rep2.rds")

MAP_estimates <- tibble(
  hc = exp(MAP_fit$par[MAP_fit$par %>%
    names() %>%
    grepl("log_hc", .)]),
  hcp = MAP_fit$par[MAP_fit$par %>%
    names() %>%
    grepl("hcp", .)]
) %>%
  mutate(municipality_idx = seq_along(hc)) %>%
  mutate(Municipality_code = municipality_idx %>% fit_data$Municipality_code_converter())

res <- readRDS("Results/vb_estimated_age_sex_structure_shrunk_hazard_rate_2021_01_18_rep1.rds")
vb_fit_age_sex <- res$vb_fit
data_age_sex <- res$data
MAP_fit_age_sex <- readRDS("Results/MAP_estimated_age_sex_structure_shrunk_hazard_rate_2021_01_18_rep1.rds")

MAP_estimates_age_sex <- tibble(
  h0m = exp(MAP_fit_age_sex$par[MAP_fit_age_sex$par %>%
    names() %>%
    grepl("log_h0m\\[", .)]),
  h0m_male = h0m * exp(MAP_fit_age_sex$par["beta_sex"]),
  nu_m = MAP_fit_age_sex$par[MAP_fit_age_sex$par %>%
    names() %>%
    grepl("nu_m", .)]
) %>%
  mutate(municipality_idx = seq_along(h0m)) %>%
  mutate(Municipality_code = municipality_idx %>% data_age_sex$Municipality_code_converter()) %>%
  mutate(Département = str_sub(string = Municipality_code, start = 1, end = 2)) 

MAP_estimates %>%
  mutate(excess_hazard = hcp, baseline_mortality_both = hc) %>%
  select(Municipality_code, excess_hazard, baseline_mortality_both) %>%
  left_join(MAP_estimates_age_sex %>%
    mutate(
      prevalence = nu_m,
      baseline_mortality_female = h0m,
      baseline_mortality_male = h0m_male
    ) %>%
    select(Municipality_code, prevalence, baseline_mortality_female, baseline_mortality_male)) %>%
  left_join(full_data %>%
    select(Municipality_code, Population, year) %>%
    arrange(year) %>%
    pivot_wider(names_from = year, values_from = Population, names_prefix = "Population")) %>%
  mutate(Département = str_sub(string = Municipality_code, start = 1, end = 2)) %>%
  write_csv("Results/MAP_estimated_hz_infect_2021_01_18.csv")


res <- readRDS("Results/vb_estimated_shrunk_hazard_rate_2021_01_18_second_round_rep1.rds")
vb_fit <- res$vb_fit
fit_data <- res$data
MAP_fit <- readRDS("Results/MAP_estimated_shrunk_hazard_rate_2021_01_18_second_round_rep1.rds")

MAP_estimates <- tibble(
  hc = exp(MAP_fit$par[MAP_fit$par %>%
                         names() %>%
                         grepl("log_hc", .)]),
  hcp = MAP_fit$par[MAP_fit$par %>%
                      names() %>%
                      grepl("hcp", .)]
) %>%
  mutate(municipality_idx = seq_along(hc)) %>%
  mutate(Municipality_code = municipality_idx %>% fit_data$Municipality_code_converter())


res <- readRDS("Results/vb_estimated_age_sex_structure_shrunk_hazard_rate_2021_01_18_second_round_rep1.rds")
vb_fit_age_sex <- res$vb_fit
data_age_sex <- res$data
MAP_fit_age_sex <- readRDS("Results/MAP_estimated_age_sex_structure_shrunk_hazard_rate_2021_01_18_second_round_rep1.rds")

MAP_estimates_age_sex <- tibble(
  h0m = exp(MAP_fit_age_sex$par[MAP_fit_age_sex$par %>%
                                  names() %>%
                                  grepl("log_h0m\\[", .)]),
  h0m_male = h0m * exp(MAP_fit_age_sex$par["beta_sex"]),
  nu_m = MAP_fit_age_sex$par[MAP_fit_age_sex$par %>%
                               names() %>%
                               grepl("nu_m", .)]
) %>%
  mutate(municipality_idx = seq_along(h0m)) %>%
  mutate(Municipality_code = municipality_idx %>% data_age_sex$Municipality_code_converter())

MAP_estimates %>%
  mutate(excess_hazard = hcp, baseline_mortality_both = hc) %>%
  select(Municipality_code, excess_hazard, baseline_mortality_both) %>%
  left_join(MAP_estimates_age_sex %>%
              mutate(
                prevalence = nu_m,
                baseline_mortality_female = h0m,
                baseline_mortality_male = h0m_male
              ) %>%
              select(Municipality_code, prevalence, baseline_mortality_female, baseline_mortality_male)) %>%
  left_join(full_data %>%
              select(Municipality_code, Population, year) %>%
              arrange(year) %>%
              pivot_wider(names_from = year, values_from = Population, names_prefix = "Population")) %>%
  mutate(Département = str_sub(string = Municipality_code, start = 1, end = 2)) %>%
  write_csv("Results/MAP_estimated_hz_infect_2021_01_18_second_round.csv")